arXiv: 1503.00655v3 [math.NA] 21 Mar 2015 


Approximation of the Scattering Amplitude using 
Nonsymmetric Saddle Point Matrices 

Amber S. Robertson a , James V. Lambers a 

a Department of Mathematics, University of Southern Mississippi, 118 College Dr #5045, 

Hattiesburg, MS 39406, USA 


Abstract 

In this paper we examine iterative methods for solving the forward (Ax = b) 
and adjoint {A 1 y = g) systems of linear equations used to approximate the 
scattering amplitude, defined by g T x = y T b. Based on an idea first proposed 
by Gene Golub, we use a conjugate gradient-like iteration for a nonsymmetric 
saddle point matrix that is constructed so as to have a real positive spec¬ 
trum. Numerical experiments show that this method is more consistent than 
known methods for computing the scattering amplitude such as GLSQR or 
QMR. We then demonstrate that when combined with known precondition¬ 
ing techniques, the proposed method exhibits more rapid convergence than 
state-of-the-art iterative methods for nonsymmetric systems. 

Keywords: nonsymmetric saddle point matrix, conjugate gradient method, 
scattering amplitude 


1. INTRODUCTION 

1.1. The Scattering Amplitude Problem 

The core objective of this paper is to design and implement an iterative 
method for the solution of a system where the coefficient matrix is large, 
sparse, and nonsymmetric. The proposed method should be more efficient 
and robust than existing methods for solving such systems. One application 
in which such a system arises is in the computation of the scattering ampli¬ 
tude. The scattering amplitude, in quantum physics, is the amplitude of the 
outgoing spherical wave relative to that of the incoming plane wave [7]]. It is 
useful when it is of interest to know what is reflected when a radar wave is 
impinging on a certain object. The scattering amplitude can be computed by 
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taking the inner product of the right hand side vector g of the adjoint system 


A T y = g 


( 1 ) 


and the solution x of the forward system 


Ax = b. 


( 2 ) 


Applications of the scattering amplitude come up in nuclear physics [l], quan¬ 
tum mechanics [HI, and computational fluid dynamics (CFD) [4|. One par¬ 
ticular application is in the design of stealth planes [if. 

The scattering amplitude g 1 x = y T b creates a relationship between the 
right hand side of the adjoint system and the solution to the forward system 
in signal processing. The held x is determined from the signal b in the 
system Ax = b. Then the signal is received on an antenna characterized by 
the vector g which is the right hand side of the adjoint system A T y = g, 
and it is expressed as g 7 x [7]. We are interested in efficiently approximating 
the scattering amplitude. It is informative to look at methods that other 
researchers have used to solve this problem, which will be discussed below. 

The solution of the linear system (J2D is important for many applications 
beyond the scattering amplitude, such as in the numerical solution of PDE 
with non-self-adjoint spatial differential operators. This solution can be ob¬ 
tained in many different ways, depending on the properties of the matrix A. 
The LDL T factorization can be used to solve some problems with a symmet¬ 
ric matrix or a Cholesky factorization can be used if the matrix is also known 
to be positive definite [9j. However, for large, sparse systems, an iterative 
method is preferred. The conjugate gradient method is the preferred iterative 
method for a symmetric positive definite matrix A [9|. However it is much 
more difficult to find this solution for a matrix that is not symmetric positive 
definite. In the case that we have a matrix that is not symmetric, we can use 
methods like the biconjugate gradient (BiCG) [3| and generalized minimal 
residual (GMRES) methods [l7| . If we have a matrix that is symmetric but 


indefinite, SymmLQ |24|, |19] is the iterative method of choice. Since the scat¬ 


tering amplitude depends on both the forward and adjoint problem, it makse 
sense to use methods that take both the forward and adjoint problems into 
account, like the quasi-minimal residual (QMR) [l6[ and generalized least 
squares residual (GLSQR) methods[25|. 
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1.2. Approximation of the Scattering Amplitude 

The method of this paper employs a conjugate gradient-like approach 
since, for large, sparse matrices, it is best to use an iterative approach, such 
as the conjugate gradient method Q which is particularly effective for sym¬ 
metric positive definite matrices. In particular, conjugate gradient has a very 
rapid convergence if A is near the identity either in the sense of a low rank 
perturbation or in the sense of the norm. In |9| it is stated that 

Theorem 1. If A = I + B is an nxn symmetric positive definite matrix and 
rank(B)—r then the Hestenes-Stiefel conjugate gradient algorithm converges 
in at most r + 1 steps. 


Theorem 2. Suppose A e M nxn is symmetric positive definite and b G M. 
If the Hestenes-Stiefel algorithm produces iterates x^, and n, = k 2 (^4) then 


x-x fc || A < 2||x — x 0 ||a 


( yHj — 1 

\ \/ft + 1 


k 




where ||w||^ = Vw T Aw. 


It is also stated in Q that the accuracy of x^. is often better than this theorem 
predicts and that the conjugate gradient method converges very rapidly in 
the A-norm if k 2 (A) « 1, where ^2(^4) is the condition number of A, defined 
by 


k 2 {A) 


PI|2P- 1 || 2 


^max(4) 

C r min( y 4) 


with (j ma y and cr min referring to the largest and smallest singular values, 
respectively. 

Multiplying both sides of Ax = b by A T yields the normal equations with 
a symmetric matrix A T A that is also positive definite when A is invertible. 
However, this approach is not conducive to solving the forward and adjoint 
problems simultaneously. Furthermore, a significant problem with using A T A 
is that now the condition number in the two-norm is squared for A T A. Since 
this increases the sensitivity of the matrix, possibly making it ill-conditioned, 
this paper explores an alternative approach. The idea is to transform the 
problems Hx = b and A T y = g into an equivalent system in which the matrix 
can be guaranteed to have real, positive eigenvalues, as well as eigenvectors 
that are in some sense orthogonal, which is then conducive to solution using 


3 





a conjugate gradient-like iteration. It is not necessarily symmetry that we 
seek, but we will have symmetry with respect to some inner product. To 
this end, we use an idea first proposed by Gene Golub in 0], and consider a 
nonsymmetric saddle point matrix that has the form 


M 


' A T WA A T ' 
-A 0 


As required by the definition of a nonsymmetric saddle point matrix, we 
assume that the matrix W is symmetric positive definite. The goal is to 
choose W so that we can guarantee M has real, positive eigenvalues. In this 
paper we will introduce the nonsymmetric saddle point conjugate gradient 
(NspCG) method to solve a nonsymmetric, large, sparse linear system, which 
will then allow us to compute the scattering amplitude. We will also use ILU 
preconditioning with NspCG, which gives rapid convergence compared to 
existing methods for solving such systems. 

This paper is organized as follows. In Section [2] we discuss the known 
methods for solving a large linear system with iterative approaches to com¬ 
pute the scattering amplitude such as Bidiagonalization or least squares 
QR (LSQR), quasi minimum residual (QMR), and block generalized LSQR 
(GLSQR). In Section [3] we will introduce the method of this paper, NspCG. 
Section 4 will include an analysis of the numerical results. The precondition¬ 
ing techniques and results can be found in Section 0 The conclusions and 
discussion of possible future work will be given in Section [6j 


2. Methods for Solving the Linear Systems of the Forward and 
Adjoint Problems 


2.1. QMR approach 

The QMR approach 00 is based on the spectral decomposition A = 
XDX -1 ] also the basis of the QMR approach is the unsymmetric Lanczos 
0,0 process which generates two sequences 


V k = [ vi v 2 ... v fc ] 

W k =[ Wi w 2 ... w fc j 

that are biorthogonal, meaning V Wk = I■ We have the following relations: 


AVk — Vk+iTk+^ki (3) 

A T W k = W k+1 T k+lik . (4) 
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where the tridiagonal matrices 


Ti 


k-\-l,k 


oil 71 

Pi «2 72 

02 ■■■ 


/3fc-i 


7fc-i 

^k 

(3k 


Tk,k 


and 


T, 


k-\-l,k 


oil 71 

01 &2 72 

A 

'' • 7fc-i 
(3k -1 

fik 


Tk,k 

hel _ 


have block structures in which and T kk are not necessarily symmetric. 
The residual, r = b — Ax, in each iteration can be expressed as 


||r fc || = || b — Axfc|| 

= ||b - Ax 0 - AVfcCfcH 

111*0 Rfe+l-^fe-l-1 ,feCfc || 

= ||14 + i(||r 0 ||ei - T k+ljk c k )\\ (5) 

with a choice of Vi = y^y where r 0 = b — Ax 0 and x/, ; = x 0 + 14cfc. We 

now have the quasi-residual ||r^|| = HHrdle! — T k+ i k c k \\. Then we choose 
w i = nfrrq where s 0 = g - A T y 0 and y k = y 0 + w k d k . Then the adjoint 

||°(J || 

residual is ||s^|| = HHsoHe! — T fc+1)fc d fc ||. The vectors c k and d fc are the 
solutions of the least squares problems for minimizing ||r^|| and ||sjj?||. So 
now the solutions can be defined as 


x fc = x 0 + V k c k (6) 

y k = yo + u k d k . (7) 


2.2. LSQR approach 

In LSQR [l, fiol ] . a truncated bidiagonalization is used in order to solve the 
forward and adjoint problems approximately. The bidiagonal factorization of 
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A is given by A = UBV T where U and V are orthogonal and B is bidiagonal. 
Thus the forward and adjoint systems can be written as 


UBV t k = b 

(8) 

VB T U T y = g. 

(9) 

Now we can solve (JBj) by solving the following two systems 


Bz = U T b 

(10) 

x = V T z, 

(11) 

and we can solve dHJ) by solving 


B t w = V T g 

(12) 

y = U T w. 

(13) 

We need to use the following recurrence relations in an iterative 
produce a bidiagonal matrix 

process to 

AV k = U k +\B k 

(14) 

A T U k +i = V k Bl + a k+ \v k+ iel +1 

(15) 


where 14 and U k are matrices with orthonormal columns, and 


Bk 


a i 

Pi «2 

Ps ■■■ 

Oik 

Pk+1 . 


Also we have that 

A T AV k = A T U k+ iB k = {V k Bl + a k+1 v k+1 el +1 )B k 

= V k BlB k + a k \ k + ie£ +1 (16) 


and 

Oik+l = Ol k+ iPk+1- 

Because B k is bidiagonal, it follows that B"[B k is symmetric and tridiag¬ 
onal. It can be seen from (UBD that m and (fT3|) implicitly apply Lanczos 
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iteration to A T A. Now this iterative process can be used to obtain the 
approximate solution to the forward and adjoint systems. We define the 
residuals at step k as 

r k = b — v4x fc (17) 

Sfc = g - A T y k (18) 


where 


x fc = x 0 + V k z k y k = y 0 + U k+1 w k , 


The goal of the LSQR approach is to obtain an approximation that minimizes 
the norm of the residual. That is, the norm |||| = ||b — Ax^H is minimized. 
When working with the forward and adjoint problems, this approach is lim¬ 
ited due to the relationship between the starting vectors 


/L T U! = ttiVx. 


The above relationship does not allow Vi to be chosen independently. 


2.3. Generalized LSQR (GLSQR) 

The GSLQR method [3, 25] overcomes the disadvantages of the LSQR 
method by choosing starting vectors Ux = and Vx = -p^ independently 
where, for an initial guess of x 0 and y 0 , r 0 = b — 7Lx 0 and s 0 = g — A T y 0 . It 
is based on the factorizations 


AV k — Uk+iTk+i'k — U k T ktk + (3 k+ iu k+ ie k (19) 

A U k Vk+iSk+i,k V k S k , k T ^] k +i^ k +i^ k (20) 

From the above we get that 

(3 k+1 u k+1 = Av k - a k u k - 7 fc_iUfc_i = c k ( 21 ) 

%+iVfc+i = A T u k - S k v k - 6> fe _xv fe _i = d fc , (22) 


where the recursion coefficients a k , y k , rj k , and 6 k are chosen to make U k and 
V k have orthonormal columns, which yields 


a k = u (23) 

Ik = ul_ x Av k+1 , (24) 

5 k = wlA T u fc , (25) 

Ok = v^iifc+x. (26) 
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We can define Ufc +1 = ^ and v*, = —, where fd k = ||cJ|, and in = II 

Pk Vk 

Now we have that 


«1 71 


i- 

Q? 

@2 OL2 


V2 d 2 ' • • 

’ ’ • ’ ’ • lk-l 

Sk+l,k 

’' • ’' • 0 k - 1 

fik 


Vk d k 

Pk+1 _ 


Vk+ 1 . 


The residuals can be expressed as follows 


Ml = || r 0 - U k+1 T k+ i >fc x fc || = |||| r 0 1|e! - T k+ljk x k \\, 


and 


Il s fc|| — || s o — VkSk+ij<yk — «fc+iVfc+ie^ + 1 y fc ||. 

The solutions x*, and y k are 


x t = x 0 + ||r 0 ||V )c T 1 ,je, 

yt = yo + Ilsoll^S^e,. 


(27) 

(28) 

(29) 

(30) 


3. Iterative Methods for Nonsymmetric Saddle Point Matrices 


The matrix M, defined as follows 


M = 


A T WA 

-A 


0 ’ 


(31) 


where A e M nxn is invertible and W is a symmetric positive definite matrix, 
is an example of a nonsymmetric saddle point matrix. It can be shown that 


x T Mx > 0 for all x/ 0. To see this, we first let x 
can be written as 


y 

z 


. Then x r Mx 


x t Mx 


' A T WA 


y 

-A 

0 

z 


y T (A T WA)y - z r Ay + y r A T z 
y T (A T WA)y. 


Now, if we let r — Ay for any nonzero vector y, then, r r — (Ay) T — 
y T A T . since W is symmetric positive definite, we have that y T (A T WA)y = 
r T Wr > 0, since r is nonzero due to A being invertible. On the other hand, 
if we assume y = 0, then x T Mx = y T [A T WA)y = 0. That is, whether y is 
nonzero or not, x r Mx = r T bbr > 0. 



3.1. Ensuring a Real Positive Spectrum 

We want to choose W so that the matrix M has a real positive spectrum, 
so it is suitable for a conjugate gradient-like iteration 
choice we need to first define 


15 . To make this 


M(i) = Jp(M) = 3{M - 7 /) = 


A T WA — 7/ A T 
A 7/ 


where p is a polynomial of degree one in the form p(() = f — 7 for 7 6 M and 


3 = 


I 0 
0 -/ 


The goal here is to determine if there exists a symmetric positive definite 
matrix M.{^) with respect to which M is symmetric, meaning that M is 
M ^-symmetric if = M t A4(j) = (M.(p/)M) T . 

Let us first define a generic nonsymmetric saddle point matrix 


A 


A B T ' 
-B C 


and then define _M( 7 ) = Jp{A). We can use the following results from 
to determine how to obtain a real positive spectrum: 



Lemma 3. Let the matrix 


J = 


I 0 
0 -/ 


be conformally partitioned with A. Then 

(1) A is 3-symmetric, i.e., 3A = A33 = (3A) T , and for any polynomial 

P, 

(2) p(A) is 3-symmetric, i.e., 3p{A) =p(A I )3 = ( 3p(A)) T , and 

(3) A is 3p(A)-symmetric, i.e., ( 3p(A))A = A T (p(A) T )3 = (3p(A)A) T . 


Theorem 4. The symmetric matrix A4( 7 ) is positive definite if and only if 

^min^Al) > 7 ^ ^max(3) (32) 

where X min mid \ max denote the smallest and largest eigenvalues, respectively, 
and 

||( 7 1 ~ C)~ 1/2 B(A - 7 /)- 1/2 || 2 < I- (33) 
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A sufficient condition that makes M. ( 7 ) positive definite can be derived from 
the above theorem. 


Corollary 5. The matrix M.[^) is symmetric positive definite when K3E) 
holds, and, in addition, 


\B\\l < (A min (A) - 7 ) (7 - A max{C)). 


(34) 


For 7 = 7 = A (A mi n {A) + X m ax{C)), the right hand side of is maximal 
and T3A ) reduces to 


2\\Bh < (Amin(A.) — A ma . T (C)). 


(35) 


The preceding results lead to a simple approach to determining whether A 
is suitable for a conjugate gradient-like iteration fljj ]. 

Corollary 6. If there exists a 7 6 I so that M.{p)) is positive definite, then 
A has a nonnegative real spectrum and a complete set of eigenvectors that 
are orthonormal with respect to the inner product defined by In case 

B has full rank, the spectrum of A is real and positive. 

Using the previous results from [lij, we obtain a simple criterion for deter¬ 
mining whether the matrix M from (l3Tf can be constructed in such a way 
as to satisfy the criterion in Corollary [ 6 ] 


Theorem 7. Let A be an invertible n x n real matrix, and let W be a sym¬ 
metric positive definite n x n matrix that satisfies 

<WW0 > fAff (36 ) 

Then the matrix M defined by 

, r f A T WA A T ' 


has real positive eigenvalues and eigenvectors that are orthogonal with respect 
to the inner product defined by A4(y) = Jp(M). That is, the above selection 
of W makes the matrix M suitable for a conjugate gradient-like iteration. 
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Proof: We need to satisfy (132]) with a proper selection of 7 . Let 


7 = i(A min (A T WA)), 

based on Corollary [5] Because of how 7 is defined, 7 satisfies 

A min (A T WA ) > 7 > 0, (37) 

which means (132]) is also satisfied. Now we need to choose W so that (133]) 
from Corollary [3] holds. We require 

2||71 t || 2 < \ min (A T WA), (38) 


or 

2o‘ max (A) < A m i n (7l T iy2L) (39) 

where ||AL t ||2 = ||A ||2 is equal to the largest singular value of A, a max (A). 
From the fact that A T WA is symmetric positive definite, we obtain 


1 

A m in(A T W A) 


= p^WA)- 1 ) 

= \\(A T WA)-% 

< l | Al- 1 ||^||^” 1 || 2 

1 

< - . 

^min (A) 2 a m i n (W) 


Therefore, (]39j) is satisfied if 




(40) 


or, equivalently, if (1361) is satisfied. □ 

It follows that the matrix W satisfies the requirements to make A4 ( 7 ) be 
symmetric positive definite and that A = M has a real, positive spectrum 
from Corollary [ 6 l This result makes the matrix suitable for a conjugate 
gradient-like iteration, as will be described below. 


3.2. The Case W = wl 

Let A = UTjV t be the SVD of A, where 

U = [ ui • • • u n ] , V = [ vi • • • v„ ] 
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and £ = diag(ui,..., a n ). In the case W = wl for some scalar w, the 
condition from Theorem [7] reduces to 


2k 2 (A) 


(41) 


We now study the eigensystem of M. Let Mx,,- = A jx.j for j = 1,2,, 2n, 

T 

where Xj = [ yj zj ] . The form of M from (l31f , with W = wl, yields 


wA T Ayj + A t zj = Xjyj, (42) 

-Ayj = XjZj (43) 

for j = 1, 2 ,..., 2 n. Substituting (J43]) into (1421) yields 

(1 - wXj)A T Zj = Xjyj. (44) 

Multiplying through by A and applying (1451) . we obtain 

(wXj - l)AA T Zj = X)zj. 


It follows that each z 3 is a multiple of a left singular vector of M, and 
X 2 3 /{wX j — 1) is the square of the corresponding singular value. Furthermore, 
from (j43j) . we find that y j is a multiple of a right singular vector of M. 

We conclude that the eigenvectors xi, x 2 ,..., x 2ri of M are given by 


X 2.7-l = 


- A . + w 

^' u i 


X 2.7 = 


- A i v i 
VjUj 


j = l,2,...,n, (45) 


with corresponding eigenvalues A = Xj~, Xj that satisfy the quadratic equa¬ 
tion 

A 2 — a 2 wX + a 2 = 0. (46) 

It can be shown directly from (1451) and (1461) that these eigenvalues are real 
and positive, and the corresponding eigenvectors linearly independent, if and 
only if w satisfies the weaker condition 


w > —, (47) 

which is consistent with the necessary and sufficient condition for A4(y) to 
be positive definite given in Theorem [41 
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3.3. Nonsymmetric Saddle Point Conjugate Gradient Method 

Let A G M nxn be nonsymmetric. We will now introduce a Conjugate 
Gradient (CG) approach that solves the linear system Ax = c by solving an 
equivalent system of the form M z = b, where 


M = 


A T WA 
-A 


A T ' 
0 


(48) 


The matrix M is also not symmetric; however, the spectrum is entirely con¬ 
tained in the right half of the complex plane, due to the fact that x T Mx > 0 
for all x. In the preceding discussion, we established that if W was chosen 
so as to satisfy the assumptions of Theorem [71 then M is diagonalizable with 
real, positive eigenvalues. Furthermore, the bilinear form (u, v)g = u J Gv, 
where G = M( 7 ) = Jy{M\ is a proper inner product, as G is symmetric 
positive definite. It follows that M is G-symmetric and G-definite, meaning 
that (Mu, v)g = (u, Mvg) for all u, v G M 2n , and (u, Mu)g > 0 for all 
u^O. 

Let the vectors p and b be defined by 


' A T Wc + d ' 


" d ' 

— c 

, P = 

0 


where Ax = c, Alz = b, and p T z = d 7 x is the scattering amplitude for 
given vectors c and d that represent the held and antenna, respectively. 
The following conjugate gradient method is based on a given inner product 
(u, v)g = v T Gu for solving the linear system Mx = b. 

Algorithm 3.1 

Input: System matrix M, right hand side vector b, inner product matrix 
W, initial guess x 0 


Require: ro = b — Mxo 

for z = 0,1 ,... until convergence do 

= (x-X,.p,) G 
1 (Pi,Pi)G 
Xi+i ^7 + cpP) 

r*+i = iy - otiMpi 

o _ (r i+ i.ppG 
W " 1 (Pi,Pi)G 

Pi+i = r i+l + Pi+lPi 

end for 

We have the inner product matrix G = At ( 7 )M suggested by [l5|. From 
15(, we see that this choice of G gives a working CG from the following 


lemma. 
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Lemma 8 . Suppose that the symmetric matrix A4( 7) is positive definite. 
Then Algorithm 3.1 is well defined for M and G = , and (until 

convergence) the scalars a* and fii+i can be computed as 


(r«, rj)A4(7) 
(Mpi, p Om^) 

**i+l)At (7) 


(50) 

(51) 


With this choice of inner product matrix, it can be shown that the residuals 
computed using the preceding algorithm are, in some sense, orthogonal. 


Theorem 9. Each residual r k as defined in Algorithm 3.1 is orthogonal to 
all previous residuals with respect to ,i.e. (rf , r fi)M(n) — 0? where i j. 

Proof: We know that r i+] = r,;—ajMp.j. Let cc* be defined as in (l50lh Also, we 
know that all of the search directions are orthogonal, i.e. pjA4fy)Mpj = 0 
for i j. We want to show that rj.M( 7 )rj = 0. This will be shown by 
induction, where the base case that we need to establish is 


- 0, i — 0,1,.... 


(52) 


To show this we use the definition of a* and the expression for the search 
directions in the above algorithm, r i+1 = r, — cpAifp t . Now we have that 




= r, - p T J ^r iM( 7 ) p. P t M 


(53) 


Reinclexing the definition of the residual from the algorithm yields the fol¬ 
lowing expression for r t 

r i = Pi- fiiPi-i. 

Substituting this into (l53|i gives 


rf +1 Mfi))ri = rf 


pfM T M(j)pi 


Rearranging the last term in (154|) yields 


iplM+M^Pi-fiipUM'M^Pi) 

(54) 


pJ'_ 1 M T M('y)/3 i Pi = fiipj M{j)Mpi_! = 0 
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because Ad ( 7 ) is symmetric, and we already know that the search directions 
p, are orthogonal with respect to Ad ( 7 ). Now it is easy to see that the 
denominator in (1541) and the last factor in the numerator cancel leaving 

if + 1 Ad( 7 )]y = rf Ad ( 7)17 - rf Ad ( 7 ) 17 = 0 . 

Now we need to show that each residual is orthogonal to all previous residuals. 
We will do this by showing rfAd( 7 )rj_ rf = 0, where d > 1. Our induction 
hypothesis is r^_ 1 Ad( 7 )r i _ (i = 0. To show this, first shift the indices to get 
the expression 

r* = iy_i - ay_iA/pj_i. 

Rearranging the recurrence relation for the search directions yields 

U—d Pi —d Pi— 1 — dfii—d • 

Using this expression for 17 and 17 we obtain 

rfA/l(7)i7_ d = Tj , _ 1 M('y)Ti- d —a i -iPi_ 1 MM('y){pi-d~Pi-i-dPi-d) 

= r^ 1 Ad( 7 )i 7 _ d - a i _ 1 pJ'_ 1 M T M(^)p i - d + 

Q! i _iP^ 1 M T Ad( 7 )p i _i_ d ^_ ( j, (55) 

where 

= 0 

by the induction hypothesis. Now we are left with 

rfA/l( 7 )i 7 _ d = -a i - 1 pJ'_ 1 M T M('y)p i - d + a i -ip^ 1 M T M{'y)p i - 1 -d/3i-d = 0, 

where both terms are 0 due to the orthogonality of the search directions. □ 

4. Numerical Results 

In this section, we will analyze the results from the methods described 
in this paper. These methods include QMR from Section 2.2, GLSQR from 
Section 2.3, and NspCG from Section 3.1. We have duplicated the results 
from (?J for GLSQR and QMR and will compare them against the results for 
our NspCG method. 

We need to first define the following matrix, M is our nonsymmetric 
saddle point matrix 

_ r a t wa -a 1 

M ~ [ A T 0 

where W = wl is defined from HTTP . These examples are from [ 7 ]. 
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4-1. Example 1 

This example uses the matrix created by A=sprand(n,n,0.2)+speye(n) 
in MATLAB where n=100. This creates a random sparse n x n matrix, where 
0.2 is the density of uniformly distributed nonzero entries, and adds this to 
the identity. 


M**' \ , 

^- 

- 


1 



i 

1 

1 



1 

1 

1 


- NspCG 

1 


- GLSQR 

1 


. QMR 

1 



1 


0 50 100 150 200 

iterations 


Figure 1: Example 1 with the matrix A 

In Figure [1] we see that at the beginning of the iteration NspCG reaches 
a better approximation in fewer iterations than either QMR or GLSQR. Al¬ 
though GLSQR eventually outperforms NspCG, it takes about 120 iterations 
before it shows any sign of convergence at all. Then it converges rapidly. 

4-2. Example 2 

Example 2 uses the ORSIRR_l matrix from the Matrix Market collec¬ 
tion, which represents a linear system used in oil reservoir modeling. This 
matrix can be obtained from http://raath.nist.gov/MatrixMarket/, 

We see that NspCG starts out with the lowest error in the 2-norm of 
the residual. Also we see that in both Figure [2] and Figure [Q that NspCG 
is more consistent than either GLSQR or QMR. Although QMR actually 
outperforms GLSQR and NspCG, it takes about 400 iterations to do so. 
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Figure 2: Example 2 


4-3. Example 3 

First define the circulant matrix 

'01 

j= 0 "• . 

1 

. 1 0 . 

Now the matrix used in this example A=le-3*sprand(n,n, 0.2)+J, where 
n=100, can be constructed in Matlab. 

NspCG starts out steady and consistent again in this Figure [3] as we see 
in Figure [2] and Figure Q] Eventually, GLSQR converges, taking about 70 
iterations to do so, while QMR fails to show any sign of convergence. 

4-4- Example 4 

We need to first define 


D 1 = 

" 1000 


e k p,p 

d 2 = 

1 

2 




1000 




Q . 
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Figure 3: Example 3 with the matrix A 

where n = p+q and £ = diag(.Di, D 2 ). Now we can define A = UY,V T , where 
U and V are orthogonal matrices. For this example we use n = 100 and D\ £ 
]R90>90_ From H] we see that NspCG starts off with the best approximation, 
but only for about 15 iterations. Then it is overtaken by GLSQR. Also, we 
can see that QMR fails to converge at all. 

4-5. Example 5 

This example uses the same definition of Di, D 2 , and A from Example 
4. In this example we will let n = 100 again, and D\ e M 50,50 . Figure [5] 
shows the same trend we have been seeing, that NspCG is more consistent 
at the beginning than any other method. At about 65 iterations GLSQR 
outperforms NspCG, and QMR fails to converge again. 

4-6. Example 6 

This example uses the same definition of Di, D 2 , and A from Example 4. 
In this example we will let n = 1000 again, and D\ £ M 600,600 . From Figure O 
we see that NspCG shows the best results for the first 600 iterations. GLSQR 
takes many iterations to converge in this case, and QMR does not converge 

at all. 
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Figure 4: Example 4 with the matrix A 


5. Preconditioning 

According to [95] conjugate gradient has very rapid convergence for a sym¬ 
metric positive definite matrix A that is nearly identity. We need to apply 
preconditioning techniques to make our matrix M satisfy this criterion. The 
result will be that the original system is transformed into an equivalent sys¬ 
tem where the coefficient matrix is near identity. As we have seen previ¬ 
ously with conjugate gradient, preconditioning techniques can be generalized 
to the nonsymmetric case. The goal is to apply ILU preconditioning jiil ]. 
while taking into account the structure of the nonsymmetric saddle point 
matrix Mde finedin (13 1 1) . The matrix W in the (1,1) block is assumed to be 
a symmetric positive definite matrix; therefore it has a Cholesky factorization 
W = GG T . We can use the OR factorization 

G t A = QR 


to obtain the factorization M = LU, where 

R Q T G- 1 
0 Q T G~ l 

Let us define 

C = G t AR~ 1 « Q, 


L = 


R T 0 
-G~ T Q G~ T Q 


U = 


(56) 
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Figure 5: Example 5 with the matrix A 


where an incomplete QR factorization [ 20 ]] is computed from the sparse ma¬ 
trix G T A which gives G T A « QR . By finding 


L~ l 


RT t J) 

R~ T QG T ’ 


U~ l 


R- 1 -Rr l 
0 GQ~ T 


(57) 


it can be seen that the resulting preconditioned system matrix is given by 

C T C~ -C T C + C T Q~ 

C T C — Q T C -C T C + C T Q + Q T C 

The above matrix has the structure similar to that of M from (13T]) . therefore 
it is a nonsymmetric saddle point matrix that is near I. 



L-'MU- 1 = 


5.1. Example 1 

The following is Example 1 from the previous section with precondition¬ 
ing. 

5.2. Example 2 

The following is Example 2 from the previous section with precondition¬ 
ing. In Figure [H NspCG converges very rapidly in only 10 iterations. QMR 
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Figure 6: Example 6 with the matrix A 

takes over 200 iterations, but still doesn’t reach the level of accuracy that 
NspCG achieves. GLSQR does not converge at all. 


6 . Conclusions and Future Work 

The results from this paper show that the NspCG method is much more 
consistent and reliable than GLSQR or QMR. NspCG only takes a few itera¬ 
tions to make fairly significant progress while GLSQR takes many iterations 
in most cases, and QMR rarely makes any progress. If preconditioning is used 
with NspCG, as is usually done with a conjugate gradient method, we have 
provided evidence that it will dramatically accelerate convergence, compared 
to state-of-the-art iterative methods such as GMRES or BiCG that are typ¬ 
ically used to solve such systems. These results support our hypothesis that 
more rapid convergence can be achieved by solving a system that, while still 
nonsymmetric, shares essential properties with symmetric positive definite 
matrices and therefore is more suitable for conjugate gradient-like iteration. 

Future work will include relating the NspCG method to a quadrature rule, 
as in [hi, 10], that can be used to compute the scattering amplitude without 
explicitly solving the forward or adjoint problem. This has been done in [7S] 
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Figure 7: Example 1 with preconditioning 


with the symmetric matrix 


[ A 0 

in conjunction with block Lanczos iteration [SjJ, but our goal is to achieve 
more rapid convergence. Furthermore, because the forward system Ax = b 
is replaced with a system with twice as many unknowns and equations, it is 
essential to implement the iteration carefully so that the gain in convergence 
speed is not offset by the additional expense of each iteration. To that end, 
it is worthwhile to consider other choices for the matrix W instead of just a 
multiple of identity. 
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